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Abstract 

The static and dynamic properties of binary mixtures of hard spheres with a diameter ratio of obI(?a = 
0.1 and a mass ratio of ms/m^ = 0.001 are investigated using event driven molecular dynamics. The 
contact value of the pair correlation functions are found to compare favourably with recently proposed 
theoretical expressions. The transport coefficients of the mixture, determined from simulation, are compared 
to the predictions of revised Enskog theory, using both a third-order Sonine expansion and direct simulation 
Monte Carlo. Overall, Enskog theory provides a fairly good description of the simulation data, with the 
exception of systems at the smallest mole fraction of larger spheres (xa = 0.01) examined. A "fines 
effect" was observed at higher packing fractions, where adding smaller spheres to a system of large spheres 
decreases the viscosity of the mixture; this effect is not captured by Enskog theory. 
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L INTRODUCTION 

Excluded volume interactions between molecules play a major role in determining the structure 
and properties of most fluids and colloidal systems. The hard sphere model, which captures the 
essence of these interactions, has played a central role in our understanding of the properties of flu- 
ids, serving as a starting point of perturbation theories for the description of real fluids 1 . Recently, 
there has been interest in binary hard sphere mixtures, where the diameters of the two compo- 
nents are very different. These systems serve as models for nanoparticle suspensions and colloid- 
polymer mixtures. In these systems, an entropically driven depletion force 2 ' 3 drives the larger par- 
ticles to cluster. While there have been many studies on the structural (e.g., radial distribution func- 
tion) and thermodynamic properties (e.g., equation of state) of these mixture G 4 ' 5i6 ' 7i8i9ilQ ' ni12113114 ' 15 , 
there have been relatively few studies on their dynamical properties. 

Much of the previous simulation work for the dynamical properties of binary mixtures has fo- 
cused on tracer particle studie a 16il?i18 , the velocity auto correlation functions, or the self-diffusion 
coefficients 19,20 , as these are relatively computationally inexpensive to determine. These studies 
have revealed that the dynamics of the larger particles deviates significantly from both the theo- 
retical predictions of Brownian particles and of Enskog theory. Lue and Woodcock 810 examined 
the self-diffusion coefficients of size asymmetric binary mixtures of hard spheres. They found a 
"fines effect" at high densities, where the addition of smaller spheres enhances the mobility of the 
larger spheres. 

Significantly less data are available for other dynamical properties. Easteal and Woolf 21 have 
investigated the tracer diffusion coefficient for binary hard sphere mixtures. They observe an in- 
verse isotopic mass effect, where heavier tracer particles diffuse faster beyond a certain solvent 
density than lighter tracer particles. Due to the computational cost of simulating highly size asym- 
metric systems, past studies have focused on small size disparity and/or moderate mole fractions 
of colloidal particles. 

Erpenbeck 2 ^ 3 - 2 ^ provided the first complete transport study, comparing predictions from En- 
skog theory and molecular dynamics results for binary hard sphere mixtures approximating a 
Helium-Xenon gas mixture. The mutual diffusion, thermal diffusion, thermal conductivity and 
shear viscosity are given over a range of state points. Enskog theory was found to provide a fairly 
good description of the transport properties for the conditions studied. Yeganegi and Zolfaghari 25 
have investigated the thermal diffusion coefficient of binary hard spheres (for moderate size ratios) 
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using non-equilibrium molecular dynamics. They observe a minimum in the thermal diffusion 
with density and good agreement with Enskog theory. Recently, Bastea 26 has investigated the vis- 
cosity and thermal conductivity of highly asymmetric "soft-sphere" mixtures at very low volume 
fractions of the larger spheres. Enskog theory was only able to qualitatively describe the results in 
that study. 

In the present work, we perform event driven molecular dynamics simulations to study the 
static and transport properties of binary hard sphere mixtures with a diameter ratio of 0. 1 and a 
mass ratio of 0.001. One of the motivations of this work is to further explore the "fines effect" 
revealed in these systems in a previous study by Lue and Woodcock 10 . Another aim of this work 
is to quantitatively test the predictive ability of the revised Enskog theory 27 for these binary hard 
sphere systems over a broad range of conditions. The remainder of this paper is organized as 
follows. Details of the hard sphere mixture model and the relation of the transport coefficients to 
the microscopic dynamics of the system are discussed in Section HO The details of the molecular 
dynamics calculations and the direct simulation Monte Carlo solution of the Enskog equation 
are provided in Section Hill The simulation data for the static and the transport properties of the 
binary hard sphere mixtures are presented in Section [TV] and the results are compared against 
the predictions of the Enskog theory. Finally, the main findings of this work are summarized in 
Section M 

II. THEORETICAL BACKGROUND 

We consider systems consisting of additive hard spheres with differing diameters and masses. 
Spheres of type a have a diameter a a and a mass m a . The spheres are not permitted to overlap, 
and so the interaction potential u ab between a sphere of type a and a sphere of type b is given by 



where r is the distance between the centers of the two spheres, and a ab = (a a + o"b) /2. Due to 
the simple nature of this interaction potential, all properties of hard sphere mixtures have a trivial 
dependence on the temperature. 

One major advantage of the hard sphere model is the simplicity of its dynamics. The dynamics 
of hard sphere systems is driven by collisions between spheres. Between collisions, the spheres 
travel at constant velocity. The solution of the trajectory of the system then reduces to determining 
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the sequence of collisions between the spheres. These collisions alter the velocities of the spheres 
but conserve their energy and momentum. After a collision between a sphere % of type a and a 
sphere j of type b, the velocities of the spheres become v- and 

v * = v * " ^ ' r ^ r ^ 

where Vj and v,- are the velocities of the spheres immediately before collision, f ^ is a unit vector 
pointing from the center of sphere i to the center of sphere j, v^- = Vj — is their relative velocity, 
and fi ab = m a m b /(m a + m b ) is the reduced mass. 

A. Static properties 

The pair correlation functions give an indication of the average local environment of the parti- 
cles in a system. For hard sphere systems, the values of the pair correlation functions at contact 
9ab{&ab) Pl a y an important role. In particular, they are directly related to the collision rates between 
the spheres: 

9ab K + 6 ) = (Anp^uy 1 (2tt (3 K 6 ) 1/2 (3) 

where p b is the number density of spheres of type b, (3 — (A: B T) _1 , k B is the Boltzmann constant, T 
is the absolute temperature, and t ab is the mean time between which a sphere of type a undergoes 
collisions with a sphere of type b. The quantity t ab can be calculated from the number of a-b 
collisions N^ oir> that occur in a simulation of duration t 

, Ngt 

t «&- 9/v (coll) ^ 

where N a is the number of spheres of type a in the system. An advantage of molecular dynamics 
simulations over Monte Carlo simulations is that the contact values of the pair correlation functions 
can be directly calculated from the times t ab and does not require the extrapolation of the pair 
correlation to contact. 

The contact values of the pair correlation functions are also directly related to the equation of 
state of the hard sphere system: 

— = 1 + -5~ 2^ X * X b a ab9ab Kb) (5) 
P 6 a,b 
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where p is the system pressure, p is the total number density of spheres, x a is the mole fraction of 
spheres of type a, and the lowercase Latin indexes run over all species (i.e. A and B for a binary 
mixture) present in the system. 

Due to the fundamental importance of the contact values of the pair correlation functions for 
hard sphere systems, there have been many efforts to develop expressions to describe them 5 '^2^. 
One of the most popular is the Boublik-Mansoori-Carnahan-Starling (BMCSL) equation of 
state^yc), ^v^ich j s an interpolation between the virial and compressibility expressions of the 
Percus-Yevick theory 31 . This is given by 

BMCSL / +\ 1 , Q~gVb £f ^|gf (f s 

9ab { ab) 1 - £3 2(1 - £ 3 ) 2 a ab + 2(1 - £ 3 ) 3 al b W 



where £„, is defined by 



6i = jE^ (7) 



Note that the solid fraction occupied by the spheres is given by <f> = £ 3 . 

The BMCSL equation yields predictions that are generally in good agreement with simulation 
data for hard sphere mixtures over a broad range of diameters and compositions 4 . However, for 
highly size asymmetric binary systems at small mole fractions of the larger spheres (often referred 
to as the colloidal limit), the BMCSL significantly underpredicts the contact value of the pair 
correlation function between the larger spheres, as compared to simulation results^^ 2 -. 

Recently, there have been several efforts to correct this. Viduna and Smith 33 — have suggested 
a new expression, based on an empirical equation of state 

«* ("J - — 3 + 2(1 _ 6)1 6— + 6(1 _ 6) , P6 + fA)-^- (8) 

This compact expression appears to compare well with simulation results. In the case of binary 
hard sphere mixtures, Henderson et al.— have suggested further modifications to the BMCSL and 
VS equations so that the contact value of the pair correlation function between the larger spheres 
yield the correct limiting behavior as the diameters of the larger spheres become infinite 32 . Their 
expressions for the pair correlation functions (which we denote as HC2) are given by 

9 H b C b 2 K.) = 9% CSL (4b) or gl% {* + BB ) (9) 

^2„2 i _ n2 t3 3 i _ #3 

_HG2 _BMCSL /-_+ ~\ i ^2°BB 1 BB 1 / lm 

^ " 9AB ^ + (l-e 3 ) 3 (l+i?) 2 (l-6f (1+i?) 3 

(4) = <?Xi (4b) + e* - 1 - x - x 2 /2 (11) 
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where A refers to the larger spheres, B refers to the smaller spheres, R = as/® a is the diameter 
ratio, and x = 3 (&(Taa ~ 6) A 



B. Calculation of transport coefficients 

In the continuum description of fluids 35 , balance equations are typically used to relate the con- 
served properties of the system (e.g., energy, momentum, and mass) to their fluxes. To close these 
equations, constitutive relations are required. These relations link the diffusive fluxes to gradients 
in the thermodynamic properties of the system. Transport coefficients are defined through the as- 
sumption that the diffusive fluxes depend linearly on the thermodynamic driving forces, which are 
gradients of local thermodynamic properties of the system. 

There are several possible choices 35 for the thermodynamic forces X and the diffusive fluxes 
J. For NVE molecular dynamics simulations, the most convenient 22 choice is the "mainstream" 
(or "unprimed"- 2 ^ 3 ^) definition of the fluxes. These are defined as 

X a = —TV (^) X A = -^VT (12a) 

J a = AiaXa + Lab'X.b J\ = L\\X\ + Lxa^a (12b) 

b a 

where [i a is the chemical potential, and J a is the diffusive flux of species a, J a is the energy flux, 
L\\ is the thermal conductivity, L a t is the mutual diffusion coefficient, and L a \ is the thermal 
diffusivity. The transport coefficients are defined through Eqs. (fT2j) . 

The relationship between stress tensor r and the strain rate in the fluid is defined in the standard 
manner: 



j)l + ( ^7? - K ) ( V • 11 ) 1 - // 



Vu+(Vuf 



(13) 



where r] is the shear viscosity, n is the bulk viscosity, and u is the streamline velocity of the fluid. 
The quantity 1 represents the unit matrix, and the superscript T indicates the transpose of a matrix. 

The Onsager reciprocity relations (L ab = L ba and L aX = L\ a ), combined with the requirement 
that J2 a J a = (due to the definition of the diffusive flux) which implies L aa = — Ylib^a Lab, 
reduce the number of independent transport coefficients to L X \, Lax, L A a, i], and k. In the fol- 
lowing section, we discuss how these transport coefficients can be determined from equilibrium 
molecular dynamics simulations. 
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C. Einstein forms of the Green-Kubo relations 



The Green-Kubo formulas relate the time correlation functions of the microscopic fluxes di- 
rectly to the transport coefficients 1 . However, the Green-Kubo relations are an unpopular method 
for obtaining the transport coefficients from molecular dynamics simulations, as they require long 
simulation times to obtain good statistics. This is not a significant issue in hard sphere systems, 
as long simulation times are more easily accessible. For systems with particles interacting with 
discontinuous potentials, the Einstein form of the Green-Kubo relations must be used, due to the 
impulsive nature of the interaction potential. The full derivation of the these formulas are already 
available 122 , and, therefore, only the final expressions are presented here for completeness. 

The Einstein relations have the general form 

m = ^- t (w^(t)w f2 (t)) (H) 

where ip(t) is a time dependent transport coefficient, V is the volume of the system, and Wyi and 
W^2 are displacement functions corresponding to time integrals of the microscopic fluxes. The 
displacement functions for a system with zero total momentum in the microcanonical ensemble 
are given in Table HI The pair of displacement functions that correspond to each of the transport 
coefficients are summarized in Table lU In hydrodynamic regime, the transport coefficients are 
given by the infinite time limit of Eq. (fl4)) 

if) = lim if;(t) (15) 

A sample of reduced correlators for a single molecular dynamics simulation run is plotted in 
Fig.CQ The function tip(t) typically displays transient behavior for short times before changing to 
the linear, long-time regime. All the transport properties, with the exception of the bulk viscosity, 
rapidly transition to the linear regime within a few mean free times. The bulk viscosity, however, 
only slowly approaches the linear regime, and, consequently, the limiting values are difficult to 
extract. As a result, we do not present data for the bulk viscosity. 

A time correlation function of a finite sized simulation is only representative of a bulk system 
for a limited duration. Beyond the time a sound wave takes to traverse the simulation box, the 
system size begins to affect the correlation function. The sound wave traversal time is determined 
directly from the speed of sound, c. For a hard sphere system the speed of sound is given by 

c 2 = m _1 fcsT 



2Z 2 dpZ 
~3~ + l)i 



(16) 
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where Z = (3p/ p is the compressibility factor, and m = J2 a x a m a is the mean particle mass. The 
HC2 equation of state (see Eqs. ©, ©, (flOl) . and CCD) is used to estimate the speed of sound, via 
Eq. (TT6T ). Data for the time correlation functions are only collected for a duration of time shorter 
than the sound wave traversal time. 

D. Enskog theory predictions for the transport coefficients 

Revised Enskog theory (RET ) 27i36i3?i38 is an extension of the highly successful Enskog theory to 
mixtures. This is the most widely applied kinetic theory of moderately dense fluids. In the Enskog 
approximation, all pre-collision correlations between particles are ignored, save for a single static 
structural correlation function. In a homogeneous system, this reduces to the values of the various 
pair correlation functions at contact, which govern the collision rates. Given these as input, Enskog 
theory yields predictions for the transport properties through the Chapman-Enskog expansion 39 . 

The standard method to solve to the Enskog equation is to expand the one-particle distribution 
function in a series of Sonine polynomials. Erpenbeck— has compiled the (corrected) Enskog 
expressions for all transport properties, excluding the bulk viscosity, of hard sphere mixtures. 
These expressions have been combined with the table of integrals given by Ferziger and Kaper 40 
and a linear equation solver to evaluate Enskog theory to the third order in the Sonine expansion. 
We present results calculated from the BMCSL and HC2 equations to determine the effect of 
improved values for g a b{&ab) on tne predictions of the transport properties. 

E. DSMC solution of the Enskog equation 

Another method for obtaining solutions to the Enskog equations is through the use of the direct 
simulation Monte Carlo (DSMC) method. This technique was originally developed for the Boltz- 
mann equation but has recently been extended to the Enskog equatio n 41 - 42 ' 43 . In this work, DSMC 
of the Enskog equation, in the style of Bird's NTC method 44 , is used to provide results. In this 
approach, the velocity distribution of each species is approximated using a set of samples 

/ a (v,t)=AT 1 ^ ( 5(v-v J (t)) (17) 

i=i 

where J\f a is the number of samples of the velocity distribution of species a. For simplicity, in 

the following expressions we assume each sample represents a single sphere. Other choices are 
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possible; however, the difference merely affects the relative sample collision testing rates and time 
scale of the simulation. 

The probability that a sample i of species a undergoes a collision event with species b after a 
time step At ab is^ 

u ih = 4:g ab (cr+) irpbcrlb [Vij ■ kj (v y ■ kj At ab (18) 

where j is a randomly chosen sample from species b, k is a randomly chosen relative orientation 
between the samples on collision, Vy = Vj — Vj is the relative velocity, and 6 is the Heaviside 
step function. The time step At ab describes the rate at which samples in species a are tested for 
collisions with a sample of species b. For a DSMC calculation of a binary mixture, there are four 
rates, one for each pairing of the species (AA, AB, BA, and BB). 

The simplest DSMC algorithm proceeds by incrementing time to the next test for collisions 
between species a and b. Each sample i of species a is tested for an event with another sample 
j, which is randomly selected. A collision is executed with a probability given by Eq. (TTBl . This 
collision only affects sample i and not the collision partner j. This method is simple but inef- 
ficient because properties that are conserved on collision (e.g., momentum and energy) are only 
conserved on average. In addition, all samples in species a are tested at each time step, which is 
computationally expensive, even though At ab is selected to yield only a few events per time step. 

An improved algorithm, based on Bird's NTC method, executes symmetric species-species col- 
lision events simultaneously, and, therefore, there are three independent test rates for the binary 
system (Ai^, AtAB = AtsA, and At BB ). For a given time step, we assume there are a maximum 
of N^ irs = NapJ^ = AWj;™ a ^ events that may occur for each species; the quantity oj\™ ax ^ is 
the maximum observed value of cu ba , which is updated, if required, during the course of a simula- 
tion. N^ irs pairs of a and b samples are randomly selected at each time step. The probability of 
collision is then scaled to 

\{2 - 5 ab )u ah ^- (19) 

pairs 

where S ab is the Kronecker delta. If the collision is accepted, then the velocities of both samples 
are updated according to the collision rule (see Eq. ©). This conserves energy and momentum 
at all times and greatly improves the statistics of the simulation. Like Enskog theory, the DSMC 
calculations require g ab (cr^ b ) as input, however, DSMC requires no polynomial expansion to make 
the problem tractable. 
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The transport coefficients are obtained through the use of the appropriate time correlation func- 
tions, as in the full molecular dynamics simulations (see Section llTBT) . DSMC provides an attrac- 
tive method of numerically solving a kinetic equation, especially as computing power increases. 
Its results are still, however, limited by the approximations of the underlying kinetic equation. 

III. SIMULATION DETAILS 

In this work, we examine the static and transport properties of highly asymmetric binary hard 
sphere mixtures. The larger A spheres have a diameter oa and mass m^, and the smaller B spheres 
have a diameter ob and mass mg. We consider systems with 03/ &a = 0.1 and m^/m^ = 0.001, 
consistent with particles of the same density. 

Discrete potentials, such as the hard sphere model, have an important advantage over more 
complex "soft" potentials. Between collisions the spheres or molecules experience no forces and 
travel on ballistic trajectories. The dynamics can be solved analytically, and the integration of 
the equations of motion is processed as a sequence of events. Current event driven molecular 
dynamics algorithms are now quite advanced and allow the simulation of large systems for the 
long times required to extract accurate transport properties. 

A. MD Simulations 

The basic event driven algorithm used in this work to perform the molecular dynamics sim- 
ulations is fundamentally the same as the one originally described by Alder and Wainwright 46 . 
Neighbor lists and the delayed states algorithm 47 are included to optimize the calculations. These 
methods are combined with a new bounded priority queue, suggested by Paul 48 , to remove the 
system size dependence of sorting the event queue. Finally, the interactions between the largest 
spheres are removed from the neighbor list and processed separately^ to allow the use of a smaller 
cell size and reduced number of collision tests. This removal is restricted to low mole fractions of 
the larger spheres as the overhead of these removed interactions is of order 0(N 2 ) in the number 
of large spheres. 

A total of iV = 13500 spheres in a cubic box of volume V with standard periodic boundary 
conditions were used in all the simulations. The volume of the system and the relative number of 
large and small spheres (i.e., Na and Nb) were adjusted to obtain the required packing fraction 
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and composition, respectively. For each of the systems examined, the initial configurations were 
equilibrated over a period of 10 7 collisions and then run for 20 trajectories of 10 8 collisions to 
collect the collision statistics and time correlation functions. 

The time correlation functions for the various transport properties were collected over approx- 
imately 100 intervals of a mean free time, using the start time averaging method™ The last 50 
values of the correlator were fitted to a line to extract the long time limit of the transport coefficient. 

B. DSMC simulations 

DSMC simulations were performed using a total of Ma +Mb = 13500 samples of the velocity 
distribution. Each of the simulations was initially equilibrated for 10 7 collisions. The time corre- 
lation functions were then collected over 8 separate trajectories, each consisting of 10 8 collisions, 
using 100 intervals of a mean free time. The statistical uncertainty of the shorter DSMC calcu- 
lations are smaller than the uncertainties of the MD simulations because Enskog theory neglects 
dynamical correlations. 

IV. RESULTS AND DISCUSSION 

In this section, we present the results of the molecular dynamics simulations for the contact 
value of the pair correlation functions and the transport coefficients of binary hard sphere mixtures. 
A comparison of the predictions of the revised Enskog theory is also provided. All quantities are 
reported in reduced units, where the unit of mass is tua, the unit of length is a a, and the unit of 
energy is k B T. 

A. Static properties 

The variation of the pressure of the binary hard sphere mixtures with packing fraction and com- 
position is shown in Fig. |2] The symbols are the data from the molecular dynamics simulations, 
and the lines are the predictions of the BMCSL (solid) and HC2 (dotted) equations of state. These 
equations of state provide an excellent description of the simulation data, with the exception of the 
very highest packing fractions where they overpredict the pressure. These deviations, however, are 
due to the onset of freezing of the larger spheres; the single component hard sphere fluid begins to 
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freeze at a packing fraction of 0.494 50 . 

The contact values of the AA, AB, and BB pair correlation functions are plotted in Fig. [3] as a 
function of the total volume fraction of spheres for different mole fractions of the larger A spheres 
xa- The simulation results for g BB are well described by the BMCSL theory. This is in agreement 
with previous simulation studies of binary hard spheres mixtures^. The VS predictions (not 
shown) provide equally accurate predictions for gBB- 

The BMCSL predictions for qab lie above the simulation results at high density for the lowest 
mole fraction studied. The HC2 predictions are higher still, however, the error is within a few 
percent. The corrections of Henderson et al. 11 to qab are small for the systems studied. The VS 
predictions (not shown) lie between the HC2 and the BMCSL results 

For the contact value of pair correlation function between the larger spheres, the BMCSL pre- 
dictions fall significantly below the simulation results at high density for the lowest mole fraction 
studied. The HC2 predictions are exceptionally accurate, even for the smallest mole fractions of 
the larger spheres. This is due to the success of the underlying VS equation (not shown), which 
give results that are nearly indistinguishable from the HC2 equation. At (f> « 0.55, qaa {^aa) f° r 
the xa = 0.5 system decreases significantly. This also occurs in the xa = OA system at a higher 
packing fraction of <p = 0.6. It appears that the larger component has frozen while the smaller 
spheres remain fluid. 

Overall, the HC2 expression is accurate and provides excellent estimates for the contact values 
of the pair correlation functions for all the conditions studied in this work. 

B. Thermal conductivity 

The thermal conductivity of the binary hard sphere mixtures is plotted in Fig. |4Jl with respect 
to the packing fraction and in Fig. 0b with respect to the pressure. The molecular dynamics 
simulation data are given by the filled symbols. The crosses are molecular simulation data for 
single component hard spheres, taken from Ref.[sjl For single component hard sphere systems, the 
thermal conductivity increases with increasing packing fraction and pressure. The initial addition 
of smaller spheres to a system of larger spheres (i.e. decreasing xa) significantly increases the 
thermal conductivity of the mixture. At the same packing fraction, a system with a lower mole 
fraction of larger spheres will have many more particles than a system with a higher mole fraction 
of larger spheres. These additional particles enhance the ability of system to transport energy. With 
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the addition of smaller spheres to the large sphere system, we observe that the thermal conductivity 
no longer increases monotonically with the packing fraction (or the pressure). Rather, the thermal 
conductivity initially decreases with increasing packing fraction down to a minimum value, and 
then it increases. The packing fraction at the minimum increases as the fraction of smaller spheres 
increases. 

Interestingly, at packing fractions below <p ~ 0.25, the thermal conductivity of pure B spheres 
(i.e., x a = 0) is lower than the thermal conductivity for the xa = 0.01 system, while for <p > 
0.25 it is higher. This implies that at sufficiently low packing fraction (or pressure) the thermal 
conductivity of the system must have a maximum with respect to xa- Physically, this would 
correspond to a situation where the addition of larger spheres to a fluid of smaller hard spheres 
would enhance its thermal conductivity. 

The solid lines in Fig. |4] are the predictions of Enskog theory within the third order Sonine 
approximation with the BMCSL expressions for the collision rates, while the dotted lines are 
the third order Enskog predictions with the HC2 expressions. The difference between using the 
BMCSL and HC2 expressions in Enskog theory is negligible, as the collisional contribution to the 
thermal conductivity is dominated by the BB and BA interactions (see Figs.[3h and c). The open 
symbols in Fig. |4] are from DSMC calculations using the HC2 expressions for the collision rates. 
These results are nearly identical to the third order Sonine approximation, indicating the accuracy 
of the approximation and validating the DSMC code. 

The simulation results are well described by Enskog theory for the pure hard sphere systems 
(i.e. xa = and 1), as well as for mixtures with relatively high mole fractions of the larger 
spheres (xa > 0.05). At high packing fractions, the Enskog predictions deviate slightly for the 
case xa = 0.5; however, this occurs at the conditions where component A appears to freeze (see 
Fig. [3k), and the BMCSL and HC2 expressions for g a b (c^) are not applicable for solid phases. 

For xa = 0.01, Enskog theory significantly underpredicts the thermal conductivity of the sys- 
tem. This deviation may be related to the enhanced mobility of the system due to the fines effect 10 
and is a result of a dynamic process not captured by Enskog theory. Note, however, that En- 
skog theory provides good predictions for the thermal conductivity of one component hard sphere 
systems 51 ', so one expects that for vanishing amounts of the larger spheres (i.e. the limit where 
xa 0), Enskog theory should again provide a fairly good description of the simulation data. 
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C. Shear viscosity 

The shear viscosity is plotted in Fig. \5\ The viscosity of all the mixtures increases monotoni- 
cally with the packing fraction of the spheres and the pressure of the system (see Fig.[5^-c). Unlike 
for the thermal conductivity, the Enskog theory predictions for the shear viscosity using the HC2 
expression for the collision rates noticeably differ from the BMCSL results (see Fig. [5k); how- 
ever, this only occurs in regions where Enskog theory poorly describes the simulation results (see 
Fig. O* and c). Enskog theory captures the low density behavior of the viscosity quite well. For 
single component hard sphere systems, Enskog theory is known to underpredict the viscosity at 
high densities 52 , due to its inability to account for correlated collisions resulting from the caging 
of spheres at these conditions. For the binary hard sphere mixtures that we study here, the Enskog 
theory underpredicts the viscosity, in general. However, the case xa = 0.01 is an exception, where 
Enskog theory actually overpredicts the viscosity at high packing fractions. 

An interesting "fines" effect occurs in the viscosity of these systems. At low overall packing 
fractions (or pressures), the addition of smaller spheres to a system of larger spheres (i.e. decreas- 
ing x A ) increases the viscosity of the system. However, above a packing fraction of about = 0.4, 
the curves for the viscosity crossover, and the addition of smaller spheres to a system of larger 
spheres decreases the viscosity of the system. This is highlighted in Fig. [5ji where the viscosity 
is almost independent of composition at a packing fraction of = 0.4. The "fines" effect is not 
captured by Enskog theory, which indicates its origin is in dynamical correlations between parti- 
cles. In these systems, the presence of the smaller spheres leads to an attractive depletion forced 
between the larger spheres, which is entropically driven. This force may disrupt the caging of 
larger spheres 10 by forcing them into closer contact, thereby creating a more open network and 
increasing the mobility of both species. 

D. Thermal diffusion coefficient 

Figure [6]presents the thermal diffusivity of the larger spheres over a range of packing fractions 
and pressures. Because Lax is negative, the larger species tends to move towards regions of higher 
temperature. Increasing the packing fraction, the pressure, or the fraction of larger spheres in the 
system decreases the magnitude of the thermal diffusivity. This general trend is in agreement with 
previous NEMD simulations 25 . 
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The use of the HC2 expressions with Enskog theory offers no significant improvement on the 
BMCSL predictions, again due to the dominance of the small spheres in the energy transport. En- 
skog theory is in quantitative agreement with the simulation data over a broad range of conditions 
examined in this work. However, the main exception is for the composition xa = 0.01, where it 
substantially underpredicts the L A \ at the higher packing fractions. 

E. Mutual diffusion coefficient 

The mutual diffusion coefficient of the binary hard sphere mixtures is plotted in Fig. [7] The 
mutual diffusion coefficient behaves similarly to the thermal diffusivity. The displacement func- 
tions required to calculated this transport coefficient contain no potential terms, and therefore, they 
do not contain a collisional component of the flux (see Tables U and [II]). Consequently, Enskog the- 
ory performs equally well with HC2 or BMCSL contact radial distribution values. Similar to the 
results for the thermal diffusivity, Enskog theory is in quantitative agreement with the simulation 
data over most of the conditions examined, with the exception of the xa = 0.01 systems, where it 
significantly underpredicts the diffusion coefficient. 

V. CONCLUSIONS 

In this work, we examined the properties of binary mixtures of hard spheres with a diameter 
ratio of a B /a A = 0.1 and a mass ratio of m B /m A = 0.001. The BMCSL equation of state is able 
to accurately describe the pressure for all the conditions that we investigated where the system 
did not freeze. However, it underpredicts the value of g A B and g A A, especially at high packing 
fractions and low mole fractions of the larger spheres. The recently developed HC2 equation, 
however, is able to quantitatively predict these quantities. 

Enskog theory provides fairly accurate predictions for the transport coefficients of the systems 
that we studied in this work. The third order Sonine approximation and the DSMC results agree 
well with one another, both validating the DSMC code and demonstrating that the third order 
solution is sufficiently accurate over the conditions studied. At low mole fractions of the larger 
hard spheres, Enskog theory fails to capture the behavior of the transport properties, especially the 
shear viscosity. This may be due to the increased correlations in the collisions between the larger 
spheres caused by the depletion forces due to the presence of the smaller spheres. 
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DSMC provides a speed benefit over traditional molecular dynamics simulations where large 
size asymmetries and low mole fractions are computationally expensive. Unfortunately, this is 
where Enskog theory begins to break down in predicting the transport properties of the fluid. 
Extension of DSMC to other kinetic theories, such as ring theory, is necessary to capture this 
behavior, however, these techniques are yet to be developed. 

We find a "fines" effect where the addition of smaller spheres to a larger hard sphere fluid 
decreases the viscosity of the system, which occurs at packing fractions greater than about 0.4. 
This effect is not captured by Enskog theory. With the addition of fines, the thermal conductivity 
of the mixture no longer monotonically increases with the packing fraction but instead initially 
decreases with increasing packing fraction to a minimum value and then increases. In addition, at 
low to moderate packing fractions, there is a region in xa where the thermal conductivity of the 
mixture is higher than thermal conductivity of either pure species. 
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TABLE I: Displacement functions for an isotropic system required to evaluate the Einstein form of the 
Green-Kubo relationships, see Eq. (fl~4l ). 





w a 




J2k a m k v k At c - c a J2k m k \r k At c 


w A 




(Ylk \m k vlw k /\t c + ImiAvfviA 


w, 


X^At c 


[12k "ifcVfcV fc At c + rmrijAvi - lpVAtcj 



The first summation runs over all time intervals between collisions At c that occur during the simulation time t. The 
indexes i and j denote the pair of spheres that undergo collision at the end of this time interval. Note that c a is the 
mass fraction of sphere of type a. 
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TABLE II: Transport coefficients and the corresponding displacement functions. The right hand columns 
indicate which rows of Table U are used. 

Lab W a , x W b ,x 

L a \ W atX W X , X 

Lxx W X , X W X , X 

V W ri,xy W v , xy 

As the system is isotropic, the transport coefficients are averaged over all components x ^ y of the displacement 
functions. 
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FIG. 1: Time dependent transport coefficients (see Eq. (fPfl)). reduced by their infinite time result, from a 
single simulation run for a binary hard sphere system with xa = 0.01 and solid fraction = 0.1. The time 
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FIG. 2: Pressure p as a function of solid fraction <p for binary hard sphere mixtures with (Tb/&a = 0.1, 
rriB/mA = 0.001, and (i) xa = 0.01 (circles), (ii) xa = 0.05 (squares), (iii) xa = 0.1 (diamonds), and 
(iv) xa = 0.5 (triangles). The filled symbols are from molecular dynamics simulations, the lines are the 
predictions of the BMCSL (solid) and HC2 (dotted) equations of state. Data points are circled where the 
system shows signs of freezing. 
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FIG. 3: Contact value of the pair correlation function g a b (cr^ b ) between the large-large (a), large-small (b), 
and small-small (c) sphere species as a function of solid fraction (p for binary hard sphere mixtures with 
pb/&a = 0.1, ms/mA = 0.001, and (i) xa = 0.01 (circles), (ii) xa = 0.05 (squares), (iii) xa = 0.1 
(diamonds), and (iv) xa = 0.5 (triangles). The solid lines are the predictions of the BMCSL equation (see 
Eq. (©), and the dotted lines are the predictions of the HC2 equation (see Eq. (fTTI)). Simulation data points 
are circled where the system shows signs of freezim* 
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FIG. 4: Thermal conductivity L\\ as a function of solid fraction <fi (a) and pressure p (b) for binary hard 
sphere mixtures with asl^A = 0.1, mfl/nii = 0.001, and (i) xa = 0.01 (circles), (ii) xa = 0.05 
(squares), (hi) xa = 0.1 (diamonds), and (iv) xa = 0.5 (triangles). The filled symbols are from molecular 
dynamics simulations, and the open symbols are the DSMC results for the Enskog theory. The crosses are 
molecular dynamics simulations for single component hard spheres, taken from Ref. 51. The lines are third 
order Enskog theory predictions using BMCSL (solid) and HC2 (dotted) values of g a b{&ab)- Simulation 
data points are circled where systems show signs of freezing. 
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FIG. 5: Shear viscosity r\ as a function of solid fraction cj> (a-b), pressure p (c), and mole fraction xa (d) 
for binary hard sphere mixtures with obI&a = 0.1 and rriB/mA = 0.001. With the exception of (d), the 
symbols indicate a mole fraction of (i) xa = 0.01 (circles), (ii) xa = 0.05 (squares), (iii) xa = 0.1 
(diamonds), and (iv) xa = 0.5 (triangles). The filled symbols are from molecular dynamics simulations, 
and the open symbols are the DSMC results for the Enskog theory. The crosses are molecular dynamics 
simulations for single component hard spheres, taken from Ref. 51. The lines are third order Enskog theory 
predictions using the BMCSL (solid) and HC2 (dotted) predictions for 5^(0"^). 
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FIG. 6: Thermal diffusivity Lax as a function of solid fraction <fi (a) and pressure p (b) for binary hard sphere 
mixtures with cfb/cfa = 0.1, ms/mA = 0.001, and (i) xa = 0.01 (circles), (ii) xa = 0.05 (squares), (iii) 
xa = 0.1 (diamonds), and (iv) xa = 0.5 (triangles). The filled symbols are from molecular dynamics 
simulations, and the open symbols are the DSMC results for the Enskog theory. The lines are third order 
Enskog theory predictions using the BMCSL (solid) and HC2 (dotted) predictions for g a b{a^). Simulation 
data points are circled where the system shows signs of freezing. 
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FIG. 7: Mutual diffusion coefficient Laa as a function of solid fraction <fi (a) and pressure p (b) for binary 
hard sphere mixtures with cfb/cta = 0.1, m^/m^ = 0.001, and (i) xa = 0.01 (circles), (ii) xa = 0.05 
(squares), (iii) xa = 0.1 (diamonds), and (iv) xa = 0.5 (triangles). The filled symbols are from molecular 
dynamics simulations, and the open symbols are the DSMC results for the Enskog theory. The lines are 
third order Enskog theory predictions using the BMCSL (solid) and HC2 (dotted) predictions for g a b{&ab)- 
Simulation data points are circled where the system shows signs of freezing. 
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